/*
    <one line to give the program's name and a brief idea of what it does.>
    Copyright (C) <year>  <name of author>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/


/**
 * \file c_HOCO.h
 * \brief  Expression HOCO
 * \author Jos de Jong
 * \date 2007-12-22
 */


#ifndef pesHOCO_H
#define pesHOCO_H
#include <msScalarFunction.h>
#include <msKineticOperator.h>

namespace atomism
{


class pesHOCO : public msScalarFunction
{
      static msRegistrar Registrar;
      static bool  isHOCORegisteredInPython;
      static msRegister* createInstance() { return new pesHOCO; }
      
protected:
  
       void registryInPython(){
             
             #if USE_PYTHON
             if( ! pesHOCO::isHOCORegisteredInPython ) {

                using namespace boost::python;
                msScalarFunction::registryInPython();

                class_< pesHOCO, bases<msScalarFunction>, boost::shared_ptr<pesHOCO> >(
                        "pesHOCO", "2D PES for HOCO coupling reaction coordinate and torsional motion" ,
                         no_init )
                .def( "New", &pesHOCO::New , "Create a new object.") 
                .staticmethod("New")
		.def( "setKineticOperator", &pesHOCO::setKineticOperator , "Set the kinetic operator. arg2: kinetic operator");
                pesHOCO::isHOCORegisteredInPython=1;
             }
             #endif
             };

       void initialize()  { 
             msScalarFunction::initialize();
	     msTreeMapper::declareChild<msKineticOperator>( KineticOperator, 
							    msKineticOperator::New( msUnitsManager::New() ), 
							    "KineticOperator");
       };
       void update(){
	     msScalarFunction::update();
       }
public:
     
       pesHOCO() : msScalarFunction(){

	     constructVar("pesHOCO","pesHOCO","PES for HOCO");
	     setConstant("a0",a0=2.5); setConstant("a1",a1=6);   setConstant("a2",a2=11);
	     setConstant("a3",a3=1);   setConstant("a4",a4=40);  setConstant("a5",a5=1.1); 
             setConstant("a6",a6=13);
       }    

       static boost::shared_ptr<pesHOCO> New(){
	    
	     boost::shared_ptr<pesHOCO> T( new pesHOCO );
	     T->initialize(); T->getYUnit().set("kcal/mol"); 
	     T->update();
	     return T;
	}

	virtual boost::shared_ptr<msTreeMapper> clone() { 

            boost::shared_ptr<pesHOCO> clone = New(); 
	    clone->msTreeMapper::operator=(*this);					      
	    return clone;
        }  

        boost::shared_ptr<msTreeMapper> setKineticOperator(boost::shared_ptr<msKineticOperator> kinop){
	
            msTreeMapper::update( KineticOperator , kinop );
	    return mySharedPtr();
	}
	
	double evaluate(){ 
	    
	    if(KineticOperator==boost::shared_ptr<msKineticOperator>())
	      BOOST_THROW_EXCEPTION(msError("The kinetic operator is not defined, the initialized function has not been called",
					    "double pesHOCO::evaluate()",getFullId())
			            );

	    msSystem& sys = *(KineticOperator->getSystem());

            KineticOperator->setDynamicDofs();
 	    sys.computeCartCoordinates();
				
	    double theta = sys.dihedral(3,2,0,1);
 	    double rHO   = sys.separation(1,3);     
            double s	 = sys.separation(0,2); 
	    double rHO_0 = sqrt( pow(   sys.separation(0,1) * sin( (180-fabs(sys.angle(1,0,2)))*csts::pi/180. ) 
                                      + sys.separation(2,3) * sin( (180-fabs(sys.angle(0,2,3)))*csts::pi/180. )     , 2 )
                               + pow(   sys.separation(0,1) * cos( (180-fabs(sys.angle(1,0,2)))*csts::pi/180. ) 
				      + sys.separation(2,3) * cos( (180-fabs(sys.angle(0,2,3)))*csts::pi/180. ) + s , 2 ));

	    double Vtot  = a0/2.*(1-cos( 2*theta*3.14/180.)) +  a1*( pow( 2.6/rHO,6) - pow( 2.6/rHO_0,6) ) 
		         + a2 * (1-a4) / (1+exp(2*a3*( s - a5 ))) 
		         + a2 * pow(1+sqrt(a4),2) / (4*pow( cosh(-a3*( s - a5 )) , 2 ))
                         + a6 ;
	    return( Vtot );
	}
		            
      	std::string getLabel(){ return ("PES for HOCO"); }        

	boost::shared_ptr<msTreeMapper> setConstant( std::string key , double value ) { 

            msScalarFunction::setConstant(key,value);

	    if( key=="a0") a0=value; if( key=="a1") a1=value; if( key=="a2") a2=value;
	    if( key=="a3") a3=value; if( key=="a4") a4=value; if( key=="a5") a5=value;
	    if( key=="a6") a6=value;
	    return mySharedPtr();
 	}

        ostream&  print(ostream& out) const {
 
            msScalarFunction::print(out);
            output::printHeader(out,"HOCO potential energy surface");
 
            out<<"Values of the parameters: "<<a0<<", "<<a1<<", "<<a2<<", "<<a3<<", "
					     <<a4<<", "<<a5<<", "<<a6<<endl;
	    return out;
 	}

 private:
	
	msChild<msKineticOperator> KineticOperator;

	double a0; double a1; double a2; 
	double a3; double a4; double a5; 
	double a6;
};

bool pesHOCO::isHOCORegisteredInPython = 0;
msRegistrar pesHOCO::Registrar("pesHOCO", pesHOCO::createInstance);

}
#endif // pesHOCO_H






